%% 模拟退火优化Simulated annealing
clear,clc
%% 参数初始化
f=@(x)1/12*sum(x.^4)-1/2*sum(x.^2)+sum(abs(x-4));% 目标函数
lb=-1e2*ones(1,6);ub=1e2*ones(1,6);%变量范围
n=length(lb);% 变量个数
%% 模拟退火计算过程
mn=10;% 进行10次，每次初始点不同
mn1=30;% 每个初始点进行30次模拟退火过程
mn2=2e3;% 模拟退火迭代最大次数
n1=4;% 每次模拟退火过程选择随机方向个数
Tm=100;% 最大温度
gama=0.99;% 退火率
result=[];% 存储mn次初始点对应最优结果集
for i=1:mn
    tic
    result1=zeros(mn1,n+1);% 存储mn1次模拟退火过程计算结果
    x0=lb+i/mn*(ub-lb);% 初始点选取
    for i1=1:mn1
        T=Tm;% 温度初始化
        x1=x0;x2=x0;% x1可接受劣解，x2不接受劣解，但迭代时以x1为准
        r2=f(x2);% 计算x2处值
        r1=f(x1);% 计算x1处值
        for i2=1:mn2
            z=-1+2*rand(n1,n);% 随机方向
            for i3=1:n1
                z1=z(i3,:);
                z1=z1/norm(z1);% 方向规范化
                xk=x1+T/Tm*z1.*(ub-lb);
                for i4=1:n
                    if(xk(i4)>ub(i4))
                        xk(i4)=x1(i4)+1/2*(ub(i4)-x1(i4));
                    elseif(xk(i4)<lb(i4))
                        xk(i4)=x1(i4)+1/2*(lb(i4)-x1(i4));
                    end
                end
                r0=f(xk);% 计算xk处值
                if(r0<=r1)
                    x1=xk;
                    r1=r0;
                    if(r0<r2)
                        x2=xk;
                        r2=r0;
                    end
                else
                    p0=exp(-1/T);% 该部分针对不同的目标函数处理方式不
                    % 同，以便x可以接受劣解,如：
                    % p0=exp(-(r0-r1)/T);
                    z2=rand;
                    if(z2>p0)
                        x1=xk;
                        r1=r0;
                    end
                end
            end
            T=gama*T;% 退火
        end
        result1(i1,:)=[r2,x2];
    end
    result=[result;result1];
    toc
end